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ABSTRACT 

• Background and Aims Zinc uptake in roots is believed to be mediated by ZIP (ZRT-, IRT- 
like Proteins) transporters. Once inside the symplast, zinc is transported to the pericycle, where it 
exits by means of HMA (Heavy Metal ATPase) transporters. The combination of symplastic 
transport and spatial separation of influx and efflux produces a pattern in which zinc accumulates in 
the pericycle. Here, mathematical modelling was employed to study the importance of ZIP 
regulation, HMA level and symplastic transport in creation of the radial pattern of zinc in primary 
roots of Arabidopsis thaliana. 

• Methods A comprehensive one-dimensional dynamical model of radial zinc transport in 
roots was developed and used to conduct simulations. The model accounts for the structure of the 
root consisting of symplast and apoplast and includes effects of water flow, diffusion, and cross- 
membrane transport via transporters. It also incorporates the radial geometry and varying porosity 
of root tissues, as well as regulation of ZIP transporters. 

• Key Results Steady state patterns were calculated for varying medium concentrations, water 
influx and HMA levels. The experimentally shown zinc gradient was reproduced very well. 
Increase of HMA or decrease in water influx lead to loss of the gradient. The dynamic behaviour for 
a change in medium concentration and water influx was also determined showing short adaptation 
times in the range of seconds to minutes. Slowing down regulation lead to oscillations in expression 
levels, suggesting the need for a fast regulation and existence of buffering agents. 

• Conclusions The model captures the experimental findings very well and confirms the 
hypothesis that low levels of HMA4 produce a radial gradient in zinc concentration. Surprisingly, 
transpiration was found to be also a key parameter. The model suggests that ZIP regulation takes 
place on a comparable time scale as symplastic transport. 

Key words: modelling, zinc uptake, ZIP, HMA, Arabidopsis thaliana, Arabidopsis halleri, 
advection, diffusion, radial transport, root, gradient, pattern 
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INTRODUCTION 



Zinc is a heavy metal and essential micronutrient for the growth of higher green plants (Sommer 
and Lipman, 1926). As part of the functional subunits or cofactor of more than 300 proteins, among 
them the zinc-finger-proteins and RNA-polymerases, it plays an important role in many cellular 
processes. There are reports of zinc ions protecting plant cells from oxidative stress mediated by 
reactive oxygen species (ROS) (Cakmak, 2000) and may act as an intracellular second messenger 
(Yamasaki et al., 2007). 

While zinc deficiency is a wide-spread problem, high doses of zinc are also disadvantageous, as 
they induce toxicity. Such elevated levels occur especially in plants growing on contaminated soils, 
e.g. in mining or industrial areas, and lead to toxicity syndroms such as reduced growth and leaf 
chlorosis in most plants (Broadley at al., 2007). Only some specialized zinc-hyperaccumulating 
species are able to tolerate high levels of zinc without any impairment (Zhao et al., 2000). In order 
to do so, they possess mechanisms for both the increased uptake of zinc from the soil and its 
sequestration and detoxification (Macnair et al., 1999). These mechanisms implicate interesting 
applications in phytoremediation or nutritional enhancement (Chaney et al, 1997; Clemens et al., 
2002) and have therefore been in the focus of scientific research. 

Zinc is normally dissolved in water as a divalent cation and in this form taken up from the soil by 
plant roots. Mobilization by secretion of chelators and acidification of the rhizosphere might be 
necessary to dissolve zinc from soil particles prior to uptake (Clemens et al., 2002). Charged zinc 
ions are unlike water virtually unable to cross cell membranes freely (Alberts et al., 2002) and 
therefore need to be transported into the root cells by specialized transporter proteins (Guerinot, 
2000; Clemens et al., 2002). These transporter proteins are tightly regulated in order to provide a 
sufficient intracellular zinc concentration without reaching toxicity. 
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Radial transport in the root 

Water and micronutrients like zinc are taken up from the soil by root cells and transported towards 
the xylem, from where they are distributed to stem and leaves (Clemens et al., 2002). On this way, 
the substances need to pass through several tissues, the epidermis, the cortex, the endodermis, and 
the pericycle (Hanikenne et al, 2008; Fig. 1). In Arabidopsis thaliana, each of these tissues 
comprises only one layer of cells (Dolan et al., 1993). The cytoplasm of adjacent cells is connected 
by cytoplasmic bridges in the cell wall, the plasmodesmata, which may be simple channels, but may 
also have complex geometries (Roberts and Oparka, 2003). Thus, the cytoplasms of all cells form a 
symplastic continuum without membrane barriers. 

Water and ions can be transported not only in this symplastic region, but also in the cell wall, the 
apoplast, which has been found to contribute significantly to root transport processes (Steudle, 
1994). The apoplastic flow, however, is stopped at the endodermis, where suberin deposited in the 
cell wall forms the Casparian strip. This strip is virtually impermeable to water and ions, although 
some findings suggest there may also be flow across this barrier (White et al., 2002; Ranathunge et 
al., 2005). Nevertheless, most water and ions need to pass the cell membrane before the Casparian 
strip to be transported further in the symplastic pathway (Yang and Jie, 2005). Because membrane 
transport is much more selective than the apoplastic flow, this barrier is thought to function as a 
mechanism to control the uptake of nutrients. Although considerable amounts of water can cross the 
membrane freely, ions are almost completely blocked and need to be transported across cell 
membranes by specialized transporter proteins (Clemens et al., 2002). 

Zinc transporters 

There are different transporter proteins involved in the transport of zinc in plants, which can be 
categorized in three families: the ZIP (ZRT-, IRT-like Proteins) family, the HMA (Heavy Metal 
ATPases) family, and the MTP (Metal Tolerance Protein) or CDF (Cation Diffusion Facilitator) 
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family. It appears that members of the ZIP family are responsible for zinc influx into the cytosol, 
HMAs efflux zinc into the apoplast (need for xylem loading), while MTPs are involved in 
sequestration of zinc to compartments such as the vacuole (Palmer and Guerinot, 2009). 

ZIP family transporters accomplish the influx of zinc into root cells. The main transporters situated 
in the cell membrane appear to be ZIP1, ZIP2, ZIP3, ZIP9, and IRT3 (Grotz et al., 1998). ZIP4 is 
also a zinc transporter, but seems to localize to the chloroplast (Guerinot, 2000). All these 
transporters are expressed at very high levels under conditions of zinc deficiency, whereas their 
expression strongly decreases within less than two hours when zinc is added to the surrounding 
medium (Talke et al, 2006). 

The exact mechanism of this regulation is mostly unknown, but recent results have shown that the 
ZIP4 gene in .4. thaliana is regulated by transcription factors of the basic-region leucine zipper 
(bZIP) family, namely bZIP19 and bZIP23 (Assuncao et al, 2010a). These factors bind to a ZDRE 
(Zinc Deficiency Response Element), which can also be found in the upstream regions of the genes 
of ZIP 1, ZIP3, ZIP4, and IRT3. 

How the bZIP19 and bZIP23 transcription factors sense internal zinc concentrations is unclear, 
since they do not appear to have a zinc binding site. It has been proposed that there are further 
players able to bind zinc and act as inhibitors of bZIP19 and bZIP23 (Assuncao et al, 2010b). It has 
been found in other regulatory networks with bZIP transcription factors that they can be regulated 
post-transcriptionally in various ways, e.g. by protein binding, phosphorylation, etc. (Schiitze et al., 
2008). Very often, these factors act as dimers (Jakoby et al., 2002). bZIP19 and bZIP23 seem to be 
partially redundant (Assuncao et al., 2010a) and it is supposed that they preferentially form 
homodimers, but may also interact to constitute heterodimers (Deppmann et al, 2006). 

The outflux of zinc from the root to the shoot mainly depends on HMA2 and HMA4 transporters, 
which are predominantly expressed in the pericycle cells adjacent to the xylem (Sinclair et al., 
2007; Hanikenne et al., 2008). Zinc hyperaccumulator species like A halleri appear to have the 
same ZIP transporters as non-hyperaccumulators, but different HMA4 genes. Moreoever, studies in 
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different plant species have shown that hyperaccumulators possess multiple copies of HMA4 in 
their genome, resulting in higher expression levels and more efficient root-to-shoot transport of zinc 
(Hanikenne et al., 2008; Lochlainn et al., 2011). 

Models of water and solute transport 

There are several approaches to model water and solute transport through root tissues. In one 
approach, tissue parts are seen as hydraulic resistors, which can be arranged in parallel and serial 
pathways treated with Ohm's and Kirchhoff s laws like resistors in an electric circuit (Steudle and 
Frensch, 1996; Steudle and Peterson, 1998). Katou and Taura (1989) and Taura et al. (1988) use 
canal advection-diffusion equations to describe water and solute movement in the apoplast. Many 
modelling approaches concern the interface between soil and root surface (Ptashnyk et al., 2011; 
Leitner et al., 2009; Zygalakis et al., 2011). To the knowledge of the authors, there has been no 
attempt to couple a structured transport model in the root tissue to a regulatory model for 
transporters. 

The regulatory mechanism of ZIP transporters has been subject of a recent model (Claus and 
Chavarria-Krauser, 2012). For the sake of simplicity, the root was seen as one compartment in that 
publication, disregarding the fact that uptake, symplastic transport and xylem loading involve 
several different cell types (Clemens et al., 2002; Hanikenne et al., 2008). Here we want to extend 
that model to consider the internal structure of root tissues in more detail. Geometry, symplastic and 
apoplastic transport will be coupled to the regulatory mechanism to understand the accumulation 
pattern of symplastic zinc and to find the prerequisites of moving zinc ions from the root surface to 
the xylem. 
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METHODS 



A modelling approach was used to conduct computer simulations of radial zinc transport in A. 
thaliana roots. The model consists of a coupled system of ordinary differential equations describing 
the regulation of ZIP transporters and 1-D partial differential equations describing the spatio- 
temporal evolution of concentration in the symplast and apoplast. Only a short description of the 
model is given below. The interested reader is referred to the supplemental file for a detailed 
derivation. All equations cited below are found in that file. 

Assumptions 

The root was simplified as a single radially symmetric cylinder and transport in the root was 
assumed to take place mainly in radial direction. This allowed to reduce the three-dimensional 
problem into coupled one-dimensional problems. The structure of the root along the radius is shown 
schematically in Fig. 1 . The root was assumed to be composed of the following cell types (from 
outside to inside): epidermis (ep), cortex (co), endodermis (en) and pericycle (pc). The cell layers 
extend from radius r x to r e . Surrounding the endodermis cells, a perfectly impermeable Casparian 
strip (cs) at position r c was assumed. 

Roots express several different influx and efflux proteins. For the sake of simplicity, we assumed 
the existence of only two types of transporters: influx (ZIP) and efflux (HMA4). Epidermis, cortex 
and endodermis cells were allowed to have ZIP transporters, while pericycle cells produced only 
HMA4 transporters. Following the results of (Talke et al., 2006), the expression oiHMA4 was 
assumed to be independent of the zinc concentration and was included into the model as a given 
amount of transporters. Both ZIP and HMA4 were assumed to be saturable and to follow Michaelis- 
Menten kinetics. The model uses no other type of signal beside the internal zinc concentration. 
Hence, coordination is achieved merely by zinc fluxes. 
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Cells have a complex internal structure with organelles, such as vacuoles, nucleus, etc. They are 
also interconnected by plasmodesmata, which represent a substantial reduction in flow cross 
section. To avoid the treatment of these internal structures, we regarded the cell content to be a sort 
of porous medium with a given volume fraction. Vacuoles were considered only by a reduction of 
flow cross section, i.e. these were not treated as separate structures and their role in sequestration 
was neglected. Cell walls were also assumed to be a porous medium of constant structure and 
porosity. We introduced a volume fraction for the symplast, which depended only on the radial 
position. This assumption is valid in view of the periodic structure of the root and orientation of cell 
layers (Fig. 1). The volume fraction of the apoplast was assumed to be constant, and based on the 
results of Kramer et al. (2007) it was set to have a value of 1/15. Fig. 2 presents at the bottom graph 
the volume fraction of the symplast used in the simulations. The values used were based on 
following estimates. The volume fraction in plasmodesmata is of the order of 0.15 (Rutschow et al., 
2011), while the vacuole was assumed to make up 0.8 of the cell volume. Other values could be 
used without problems. 

Water flow 

Within each compartment, zinc is carried along the water flow path (advection). This process 
influences the distribution of zinc and determines how fast variations in external zinc concentration 
spread in the system. To avoid a complete treatment of water fluxes in root tissues, we focused only 
on mass conservation delivering the radial flow speed by considering the effective flow cross 
section. Variation of cross section in the symplast was included by the volume fraction shown in the 
bottom graph of Fig. 2. Water fluxes in the apoplast were assumed to be small and were neglected. 
Hence, water exchange was assumed to occur only at the xylem-pericycle and epidermis-medium 
interfaces, although the apoplast is believed to contribute to the total flux (Steudle, 2000). 
Epidermal cells were assumed to take up water from the medium with a given constant flux velocity 
qo. Different transpiration rates were simulated by variation of qo. This approach is very simplistic 
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and other more sophisticated models have been proposed, for example (Katou and Furumoto, 1986; 
Katou et al., 1987; Taura et al, 1988; Katou and Taura, 1989; Murphy, 2000). Modelling water 
fluxes in plant tissues is a complex problem which deserves treatment of its own and is out of scope 
of this manuscript. 

A consequence of mass conservation, incompressibility of water, and the radially oriented flow is 
that the mass flux through two concentric surfaces must be equal. As the cross section decreases 
towards the stele, the flow speed has to increase proportionally to compensate for the smaller cross 
section (see Fig. S.3 in supplemental file). The top graph in Fig. 2 presents the flow velocity used in 
the sequel. It shows clearly how the cylindrical geometry of the root results in a general increase of 
velocity towards the stele. 

Zinc transport 

Transport of zinc was modelled by diffusion in the cell wall and diffusion-advection in the 
symplast. Note that we assumed no water fluxes in the apoplast, which means consequently that 
advection in the apoplast was neglected. Using mass conservation of zinc, three equations for the 
concentration in the apoplast and in the symplast apply [supplemental Eqs. (Slla,b)]. Solving these 
equations would deliver the time evolution of three dimensional distributions of zinc in the root 
tissue. For this purpose, a precise 3-D representation of the tissue and computationally expensive 
numerical methods would be needed. Therefore, considering the structure of the tissue shown in 
Fig. 1, we focused on the radial distribution by reducing Eqs. (SI la,b) into a system of 1-D 
equations [supplemental Eq. (SI 5)]. These equations describe the time evolution of the radial 
distribution of zinc in the apoplast and in the symplast, and were those used to conduct the 
simulations. 
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In addition to diffusion and advection, zinc fluxes through the membrane (ZIP and HMA4 
transporters) have to be considered. These fluxes are modelled as chemical reactions taking place 
on the membrane 

apoplastic Zn + ZIP <-> Zn-ZIP — > symplastic Zn + ZIP 
symplastic Zn + HMA4 <->■ Zn-HMA4 —> apoplastic Zn + HMA4 

The first reaction describes uptake by ZIP, while the second reaction describes efflux by HMA4. 
Because the amount of free transporters is finite, the mechanism must saturate. As proposed by 
Claus and Chavarria-Krauser (2012), saturation was modelled by a Michaelis-Menten saturation 
function instead of considering the full reactions given above. Influx was assumed to take place on 
the surface of epidermis, cortex and endodermis cells, while the efflux was assumed to take place 
only at pericycle cells. The mechanism was coupled to the zinc conservation equations Eq. (SI 5) 
via a transporter dependent zinc flux density [compare supplemental Eqs. (SI 2) and (SI 8)]. More 
details on the transport model and its derivation are found in the supplemental file. 

Regulation 

Regulation of ZIP has been discussed in detail in a previous paper (Claus and Chavarria-Krauser, 
2012), where we found an activator- inhibitor model with dimerization to present a likely 
mechanism. This model was used here to describe the amount of ZIP transporters for each cell type 
shown in Fig. 1 besides the pericycle. A scheme of the mechanism is presented in Fig. 3. An 
inhibitor / senses the internal zinc concentration and inhibits a dimerising activator A. Gene 
expression for the production of transporters depends on activation by A. An increase in internal 
zinc results in higher levels of inhibitor, lower levels of activator, lower production of transporter 
and, thus, less zinc influx - and vice versa for a decrease in concentration. More details on the 
regulation model are given in the supplemental file. 

Although this regulatory mechanism was developed as an average description over the whole root, 
it is more realistic to assume individual gene expression and transporter production in each single 
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cell. These regulatory systems had to be coupled to transport in the apoplast and symplast [Eq. 
(SI 5)], which deserves a few considerations. First, zinc is not homogeneously distributed within 
one cell and it is unclear where and how the cell actually senses the concentration. Here, we 
assumed that a cell senses the average concentration inside it. Depending on the sensed 
concentration, each cell adjusts its own expression level independently of the other cells. The 
resulting amounts of transporters for endodermis, cortex and epidermis cells need to enter the 
transport model as a spatial distribution. They are assumed to be evenly distributed on the plasma- 
membrane of each single cell, but differ from cell to cell. This assumption is supported by the 
HMA2 expression pattern found in A. thaliana (Sinclair et al., 2007) (Fig. 1(c) of that publication). 
Hence, a piecewise constant spatial distribution of transporters was coupled to the transport model 
Eq. (SI 5) via the zinc flux density given in Eq. (SI 8). The reader is referred to Eqs. (S3), (S4) and 
(S16a,b) in the supplemental file for a precise definition of the transporter distributions, while the 
sensed concentration is defined in Eq. (S2). 

Parameters 

The diffusivity of zinc ions in water has been measured by Harned and Hudson (1952) to be 700 (im 
s" 2 Values of 530 jam s" 2 have been found for calcium in the axoplasm of Myxicola (Donahue and 
Abercombie, 1987). Based on these values we assumed an approximate diffusivity of D = 500 um s" 
2 for zinc in the symplasm. This value does not consider the diffusion of chelated zinc, which could 
be several times slower because of the added mass. Due to its porosity, diffusion in the cell wall has 
been measured to be fifteen times slower than in the cytoplasm (Kramer et al, 2007), giving a 
volume fraction of 1/15. The reduction in diffusivity was accounted for correctly in the model 
through inclusion of this volume fraction [compare Eq. (SI 5)]. For an estimate of the surface water 
flux velocity q we relied on data from Rosene (1943) for water fluxes through a root hair cells and 
data from Zarebanadkouki et al. (2012) for Lupine albus roots, qo was then manually adapted to 
obtain similar patterns as observed by Sinclair et al. (2007) and Hanikenne et al. (2008). Other 
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unknown parameters, for example H and T , were also manually adjusted to qualitatively 
reproduce those data. The parameters of the regulation model were adopted from Claus and 
Chavarria-Krauser (2012), where an automatic fitting procedure was used to obtain them from 
measured expression levels. Sizes of cells and tissue layers were estimated from Hanikenne et al. 
(2008) with r e = 40 |am, r c = 12.5 um, and r x = 6 |j,m. The most important model parameters are 
listed in Table 1 and the cell sized used are listed in Table 2. Other parameters are found in Tables 
S.l to S.3 in the supplemental file. 

Numerical methods 

The system of one-dimensional partial differential equations [Eq. (SI 5)] was solved numerically 
with a conservative finite difference method. Operator splitting was used to employ different stable 
explicit finite difference schemes to the advective and diffusive contributions. For the advection 
term, a second-order McCormack method was applied. A first-order FTCS (Forward Time Centered 
Space) scheme was used for the diffusion and reaction terms. The boundary conditions were 
implemented with upwinding. To guarantee mass conservation, the boundary flux was corrected for 
numerical diffusion of the scheme. The regulation model Eq. (SI) was solved by an explicit Euler 
method. Using explicit schemes for both systems allowed to couple the solvers without much effort. 
Fig. S.2 in the supplemental file illustrates the algorithm used to couple the different numerical 
schemes. Steady states were calculated with Newton's method using numerical derivatives. Average 
gradients from spatial data were obtained by linear regression. 
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RESULTS 

Steady state 

Figure 4 shows the simulated steady state patterns of ZIP levels (boxes), and the zinc concentrations 
in the symplast (solid line) and apoplast (dashed line) for roots grown in media with high (Z e = 25 
uM) and low zinc (Z e = 1 uM). In both cases, the apoplastic zinc concentration decreased from the 
epidermis towards the casparian strip, as a consequence of uptake by ZIP, although this effect was 
clearer for high zinc. Behind the casparian strip, the apoplastic concentration rose again due to 
HMA4 mediated efflux from the symplasm. The apoplastic concentration inside the stele was 
smaller for low external zinc. However, it was still about 60% of the concentration at high zinc, 
although the external concentration was 25 times smaller. The concentration in the xylem was very 
low in both cases, because inflowing water diluted the solution strongly. Regarding the symplastic 
concentration, a continuous gradient towards the stele was found. This gradient was substantially 
less obvious at low zinc, while strongly pronounced at high zinc. ZIP transporters were almost 
completely down-regulated in all cells for high zinc, while their expression was much higher for 
low zinc. The model predicts a gradient in ZIP activity with high activity in the epidermis and low 
activity in the endodermis. Although this gradient was still present for high zinc, the difference in 
expression level was very small and is probably hard to detect in a measurement. 

Variation of HMA4 and water influx velocity 

To model A. thaliana HMA4 mutants, simulations were performed with higher levels of HMA4 
transporters (via parameter H ; Table 1). Figure 5 presents in the upper panel the symplastic 
concentration distribution for different HMA4 levels. Increasing HMA4 to two (dashed line) or 
three times (dotted line) of the original level (solid line) lead to a decrease of the overall zinc 
concentration and loss of the gradient. Surprisingly, a very similar behaviour was seen when 
varying the influx velocity of water (parameter q ; Table 1). Simulations with half (dashed line) or 
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one fourth (dotted line) of the original velocity (see Figure 5, lower panel) showed a loss of 
symplastic zinc gradient similar to increasing HMA4. The overall zinc concentration, however, 
remained higher than in the variation of HMA4, with even slightly higher values in the epidermis 
than for the original velocity (q =4 urn s" 1 ). Very low water fluxes and very high HMA4 levels 
produced even a retrograde gradient (data not shown). 

To investigate further the relation between HMA4 level and symplastic gradient, the correlation 
between steady state average gradient and HMA4 activity was analysed. The average gradient was 
calculated by linear regression as the slope of a fitted linear function. In the case of an accumulation 
in the pericycle the result is a negative number, reflecting the fact that the gradient points towards 
lower radii. The absolute value of this gradient gives a measure of the "steepness" of the 
concentration gradient. Figure 6 shows the correlation between the average gradient and the 
normalized HMA4 activity (solid line) or the external zinc concentration (dashed line). While the 
response to variation of external concentration was fairly linear with an increase in steepness 
following an increase in concentration, the response to HMA4 was strongly non-linear. Raising 
HMA4 lead to a strong decrease in gradient steepness. At three times the original level, the gradient 
became even retrograde (change in sign). At less than half the original value, a singularity/pole was 
found. The gradient became very steep when approaching the pole and no steady state can be 
sustained beyond. The reason is an imbalance between influx and efflux, leading to more zinc being 
pumped into the symplast than can efflux into the xylem. Small variation of HMA4 produced a 
strong reaction near the pole, i.e. zinc accumulation reacted very sensitive to HMA4 activity there. 

Dynamics 

The time evolution of the response of the system to changes in environmental conditions was 

simulated to understand the transient behaviour of the interplay between regulation and transport. 

Starting from the steady state at low zinc (Z e = 1 uM), the external zinc level was raised to a high 

zinc condition (Z e = 10 |aM). This procedure corresponds to a resupply experiment. Figure 7 
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presents the evolution of the distributions in time during this adaptation. The apoplast outside the 
Casparian strip adapted to the new condition very quickly, while the rest of the system reacted more 
slowly. After increasing the external concentration zinc diffused quickly into the apoplast outside 
the Casparian band (in less than 1 s). During this time, regulation kept the high expression level of 
ZIP, resulting in an overall increase in symplastic concentration and a more pronounced gradient. 
This lead to an increase in apoplastic concentration inside the stele and later to a sudden down- 
regulation of ZIP. Although the time scale of regulation was assumed to be equal to that of transport 
(parameter t = 1; Table 1), adaptation of ZIP activity lacked behind, leading to an "overshoot" at 
about 5 s, where symplastic zinc in the pericycle exceeded its final steady state value by a factor of 
two. From there on the system stabilized to reach the new steady state value after about 20 s. 

Transpiration controls the velocity of water flow, can vary substantially in the course of the day and 
is minimal during the night. To understand how adaptation to a change in transpiration rate takes 
place, we simulated the time evolution upon a sudden change in water influx velocity q from 0.05 
to 4 um s" 1 (see Figure 8). Due to the low water flux, a slightly retrograde concentration gradient 
was found at s. ZIP expression was also almost constant in space, with a slightly higher 
expression in the endodermis. After increasing q , the concentration gradient built up quickly. An 
accumulation in the pericylce was clearly visible after 0.2 s and stabilized already after 1 s. This 
lead to a reversal of the distribution of gene expression, with a clear gradient towards the epidermis, 
where expression stabilized at twice the level of endodermis cells. The apoplastic concentration 
outside the Casparian strip did not change much during equilibration. In contrast, the apoplastic 
concentration in the stele decreased to half the initial value as a consequence of the higher water 
flow rate. Also noteworthy is that due to the same reason, the concentration in the xylem fell from 
about 10 uM for low transpiration to almost uM for a normal transpiration rate. 
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Time scale of regulation 

In the experiments described above, the time scale of regulation was set to be comparable to the 
time scale of transport and diffusion (parameter t = 1; Table 1). The real time scale of regulation is 
unknown and may indeed be much slower. To understand how the specific choice of x influences 
the entire process, simulations with different time scales were performed. Figure 9 shows the effect 
of resupply from low zinc (Ze = 1 |iM; upper black square) to high zinc (Ze =10 lower black 
square) by plotting: the evolution in time of average ZIP activity and average internal zinc 
concentration as a phase diagram (Fig. 9A), and the average internal concentration against time 
(Fig. 9B). The paths shown in Fig. 9A represent the state of the root in time and the arrows mark the 
direction in which the state moved. The same transition from low to high zinc was conducted for 
three different time scale factors x: 1 (red path), 0.1 (blue path) and 0.01 (green path). This means 
that the time scale of regulation was approximately 10 s for x = 1, 100 s for x = 0.1 or 1000 s for x = 
0.01. 

In general, the internal concentration rose strongly after resupply, exceeding the one of the new 
steady state at high zinc ("overshoot", Fig. 9A, B). The reason was the high ZIP expression level at 
the initial state (about 30% activity; Fig. 9A). Regulation then reacted by shutting down the 
expression of ZIP (vertical portion of the paths; Fig. 9A), overreacting even slightly. While for x = 1 
the system eventually reached the new steady state, the cases with slow regulation (x = 0.1 and 
0.01) oscillated around the steady state on a limit cycle without reaching it (Fig. 9A, B). The reason 
for the oscillation was the overreaction of regulation producing a strong sudden reduction in 
concentration (horizontal path; Fig. 9A), which lead again to an overreaction of up-regulation, and 
to a too large concentration, etc. The amplitude of the oscillation correlated with the time scale of 
the regulation: slow regulation produced larger oscillations (Fig. 9A, B). 
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DISCUSSION 

Steady state 

Fluorescence imaging data (Sinclair et al., 2007; Hanikenne et al, 2008) show higher zinc 
concentrations in the apoplast than in the symplast. A similar distribution is seen in our simulations, 
although the accumulation in the apoplast may not appear as prominent as in the fluorescence 
images. This may be due to the fact that zinc and other cations can bind to components of the cell 
wall and accumulate in the apoplast (Sattelmacher, 2001), which is not considered in our model. 
Also, the fluorophore Zinpyr-1 used by Sinclair et al. (2007) and Hanikenne et al. (2008) reflects 
only levels of non-chelated zinc (Sinclair et al, 2007). Considering, in addition, that the vacuole 
contributes up to 90% of the cell volume, the fluorescence images reflect rather the concentration in 
the vacuole. All these points could result in underestimating the symplastic zinc concentration in 
measurements. 

Within the symplast, a radial concentration gradient with accumulation in the pericycle has been 
found in experiments (Sinclair et al, 2007; Hanikenne et al., 2008). This pattern was reproduced 
very well by the model under certain conditions. We found that for existence of the pattern the 
contributions of influx, efflux, diffusion and advection have to be properly balanced. Here, influx 
into the symplast (parameter T ; Table 1) needed to be about 100 times higher than efflux into the 
apoplast (parameter H ; Table 1) to sustain the pattern in equilibrium. Although the absolute value 
of Ho and To are more or less arbitrary, their ratio should give a good estimate for a real root. Much 
less HMA4 than ZIPs was needed to obtain the pattern, which could explain why roots express so 
many different ZIPs compared to HMAs. Higher influx or higher efflux produced retrograde 
gradients, when advection was not increased correspondingly. This result is surprising in the sense 
that the pattern can be expected to vary in the course of the day, since advection varies in roots as a 
consequence of changes in transpiration rates. Diffusion is by its nature an equilibrating process, 
which seeks to even out any concentration gradient. Hence, a large diffusivity compared to 
advection would destroy the pattern. 

17 



The two crucial points for the creation of the pattern are water flow and geometry. First, the radially 
oriented advection produced by water uptake is the only process that is physically able to create 
accumulation at the pericycle. It links the spatially separated influx and efflux cells. Although these 
are also linked by diffusion, accumulation cannot be explained by that process. The importance of 
the velocity of water influx will be discussed below in more detail. Second, the cylindrical 
geometry is beneficial to creation of the pattern, because the volume contracts for smaller radii. On 
the one hand, this accelerates water on its path to the xylem, helping to create a larger accumulation 
in the pericycle. On the other hand, the concentration increases faster at small radii for the same 
flux, because there is less solvent volume. Without the geometrical effects, the pattern was far less 
pronounced (data not shown). The influx velocity of water, q , needed to be larger than 1 um s" 1 to 
produce a sufficient gradient (Fig. 5 bottom). Zarebanadkouki et al. (2012) measured velocities of 
approximately 0.2 um s" 1 in Lupine albus during the day. The measured values are lower than the 
values used here, but these are for lateral roots of L. albus with a diameter substantially larger than 
the primary root of A. thaliana. Measurement of the flux velocity in A. thaliana roots and an 
extension of the water flow model would be necessary to be able to draw more precise conclusions. 

Regarding the patterns of ZIP expression, our model predicts that ZIP activity follows conversely 
the pattern of symplastic concentration. The gradient in expression is particularly clear for low zinc 
(Z e = 1 uM), where the epidermis is predicted to have an expression level that is five times higher 
than in the endodermis. In Claus and Chavarria-Krauser (2012) we showed that the regulation 
mechanism proposed is particularly robust for certain internal concentrations (Fig. 6A of that 
publication). Gene expression can vary substantially without affecting much the internal zinc 
concentration in this "robust" range. Our model suggests that a high dynamic range in expression 
level, i.e. high and low expression in the epidermis and endodermis, respectively, can be expected 
near 1 uM external zinc concentration. This is the external concentration at which an experimental 
validation of the expression pattern should be conducted to obtain the best results. Measurements of 
Birnbaum et al. (2003) are not in concord with our results. Using growth medium with 30 uM 
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external zinc, these authors found the expression levels of ZIP2 and ZIP4 to be minimal in the 
cortex and similarly high in the epidermis and endodermis. The measurements of Birnbaum et al. 
(2003) are surprising in the sense that they are not in line with the well documented symplastic 
gradient. Their results suggest that concentration of zinc is maximal in the cortex cells, where 
expression was found to be lowest. Neither our model nor the fluorescence images of Sinclair et al. 
(2007) and Hanikenne et al. (2008) support this. As mentioned above, an experimental validation 
of the expression patterns should be conducted at an external concentration for which a high 
dynamic range in expression is expected. At 30 uM external zinc, where Birnbaum et al. (2003) 
conducted their measurements, expression can be expected to be very low in general and the 
dynamic range per se small (see for example Fig. 4 top). Small gradients in expression level are 
probably insignificant compared to the uncertainty of the measurement at that high concentration. 
Future experiments are needed to verify the gradient in ZIP activity predicted by the model. 

Variation of HMA4 and water influx velocity 

The level of HMA4 has been increased in experiments by introducing a HMA4 gene of the zinc 
hyperaccumulator A. halleri into roots of A. thaliana (Hanikenne et al., 2008). The AhHMA4 gene 
has been found to contain multiple copies of HMA4, which leads to higher expression and more 
efficient transport of zinc into the xylem (Hanikenne et al., 2008). As a result, fluorescence images 
showed a change in the distribution of zinc in the tissue such that accumulation in the pericycle and 
the radial gradient were lost. Indeed, this effect is captured by the model very well, as increasing 
HMA4 to two or three times the original level lead to a decrease of the overall zinc concentration 
and a loss of the gradient (Fig. 5 top and Fig. 6). The contrary effect has been observed in hma2, 
hma4 double mutants of A. thaliana (Sinclair et al., 2007) and in A. halleri with reduced expression 
of AhHMA4 (Hanikenne et al., 2008). The model describes this situation also very well (Fig. 6) and 
predicts a high sensitivity to variations in HMA4 in this regime. Talke et al. (2006) showed that 
HMA4 expression in A. thaliana and A. halleri roots varies substantially less in resupply and 
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oversupply experiments than the one of ZIP genes. They used this fact to propose that its expression 
does not depend much on the zinc status. In view that our model predicts a regime of large 
sensitivity, HMA4 might actually also be subjected to regulation and only small adaptation of 
expression might be enough to create sufficient effect. This would explain why HMA4 expression 
varies less. 

Surprisingly, decreasing the influx velocity had a very similar effect as increasing HMA4 (Fig. 5). 
Since plants change their respiration rates depending on the time of day, light, humidity, etc. the 
flow velocity of water in roots is highly variable. Simulations with half or one fourth of the original 
velocity showed a loss of the radial zinc gradient, similar to an increase in HMA4. The only 
difference was the overall zinc concentration, which remained slightly higher. While an increase in 
HMA4 caused enhanced efflux of zinc into the xylem and thereby an increase in the total outflux 
from the symplastic domain, a decrease in velocity only changed the distribution of zinc in the 
tissue. Instead of accumulating zinc in the pericycle, diffusion took overhand and produced an 
almost homogeneous distribution. Adaptation to a new transpiration rate is predicted to take place 
within less than a minute (Fig. 8). No published experimental data examining the relation between 
transpiration rate and zinc localization in roots is known to the authors. If the equilibration times 
predicted by the model prove to be realistic, experiments will have to be designed and conducted 
with much care to avoid artefacts - just moving the plant into dark to conduct the measurement 
might change the pattern. 

External zinc concentrations are predicted by the model to influence the strength of the zinc 
concentration gradient in an almost linear manner (Fig. 6), i.e. higher external concentrations lead to 
higher accumulations of zinc in the pericycle. No experimental quantifications of the concentration 
gradient for varying concentrations are known to the authors. Sinclair et al. (2007) found an almost 
linear relation between Zinpyr-1 fluorescence and zinc concentration in the medium, suggesting that 
the average internal concentration depends linearly on the external concentration. Our model 
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predicts also a roughly linear dependency with a slope of 0.0887 ± 0.0005 uM uM" 1 . Talke et al. 
(2006) also measured the average zinc concentration in A. thaliana and A. halleri roots for different 
external concentration and found a positive correlation. They plotted their data with a logarithmic 
scale, which does not allow to compare directly the functional relations. However, plotting a linear 
relationship against a logarithmic scale produces a graph similar to the one published by Talke et al. 
(2006), suggesting that at least for small concentrations the relationship is roughly linear. 

Dynamics 

Similar to the resupply experiments by Talke et al. (2006), we simulated the adaptation of the 
system to a new environmental condition with increased external zinc. In our simulations, the new 
steady state was reached within 20 seconds (Fig. 7). Since adaptation of the regulatory system is 
relatively slow compared to transport, an "overshoot" was produced, where zinc concentration 
temporarily exceeded the final value before reaching the new steady state. Talke et al. (2006) 
measured the expression levels every two hours, which is too coarse to resolve the dynamics of 
adaptation in the root, as suggested by our results. 

While a sudden change in the external zinc concentration is rather unlikely in a natural 
environment, changes in the flow velocity happen regularly caused by changes in leaf respiration 
rates. Our simulations showed that here the new steady state is reached within 5 seconds and that 
there is no "overshoot" as described above. 

Equilibration predicted by our model seems at a first glance to take place very fast (within less than 
a minute). Due to the lack of data, the time scales of both transport and regulation had to be chosen 
more or less arbitrarily, which makes this result uncertain. Nevertheless, the following reflection 
indicates that equilibration within minutes is feasible. The time scale of transport is in principle 
determined by four parameters beside the diameter of the root (Table 1): q , D, H and T . While q , 
H and T were chosen such that our results correspond qualitatively to the observations of Sinclair 
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et al. (2007) and Hanikenne et al. (2008), the diffusivity D in the cytoplasm was estimated from 
measurements of free calcium in the axoplasm of Myxicola (Donahue and Abercombie, 1987). That 
is, D was set while the rest was chosen to match the observations. The time scale can be adapted by 
multiplying these parameters with the same factor. Therefore, if the diffusion coefficient D would 
be smaller, q , T and H could have been chosen correspondingly smaller obtaining the same spatial 
pattern but with a slower equilibration. Zinc is known to be chelated in the cytoplasm (Clemens et 
al., 2002). Depending on the size of the chelator, diffusion coefficients may be an order of 
magnitude lower. Hence, the time scale of our model would be roughly ten times higher and 
equilibration would take place in three to four minutes. The time scale of equilibration should be 
considered with more detail in further experiments. Our model indicates that experiments should be 
conceived to capture effects that live only a few minutes. 

Simulations with different time scales of regulation were performed to understand the origin of the 
overshoot found during equilibration. Variation of the relative time scale between transport and 
regulation showed that a slow regulation lead to oscillations with high amplitudes (Fig. 9). 
Experimental validation of these oscillations may be difficult, because a real root lacks a high 
degree of spontaneous coordination. The model assumes that the root is perfectly coordinated 
resulting in well defined oscillations. In real roots, the oscillatory behaviour would rather come to 
light as a high variance in, for example, expression level. Since zinc shocks caused by oscillations 
can be toxic and even deadly to plant cells, regulation either needs to be fast - on the same time 
scale as transport - or the cells need to have fast and efficient buffering mechanisms. So far, our 
model neglects buffering and sequestration, although zinc is known to be sequestered into the 
vacuole (Clemens et al, 2002). To be efficient, these mechanisms would have to be fast and not rely 
directly on regulation, which calls for a chelation mechanism. This point, however, needs to be 
treated in more detail in a future extension of the model. 
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Conclusions 

Water uptake - and the herewith associated advection of zinc towards the stele - is the main 
mechanism in creation of the radial zinc pattern in roots. The cylindrical geometry is a benevolent 
factor supporting the effect of advection. The transpiration rate can be expected to influence the 
pattern strongly, since a reduction should produce similar effects to increasing the expression of 
HMA4. Zinc accumulation in the pericycle depends non-linearly on the expression of HMA4, where 
small variations of expression suffice to produce large effects. This might explain why HMA4 
seems to be unchanged during resupply experiments. In general a much smaller level of HMA4 
than of ZIPs is needed to maintain sufficient zinc supply, giving a possible explanation why so 
many different ZIP transporter are expressed in roots. Resupply of zinc showed that regulation has 
to take effect within minutes to avoid strong peaks in symplastic concentration. Fast chelating 
agents seem to be necessary to dampen possible oscillatory behaviour and short-term oversupply. A 
slow reacting sequestration possibly based on genetic regulation, such as one into the vacuole, 
would be less suited to counteract this short-term effects and may be rather important in long term 
adaptation. We conclude in general that regulation has to be faster than assumed before. 
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TABLES 



Table 1 : List of parameters 



Parameter 


Value 


Description 


qo 


0.1 -4 um s" 1 


Water flux velocity at root surface 


D 


500 um 2 s" 1 


Diffusivity of zinc in cytoplasm 


r x 


6 urn 


Stele radius (without pericycle cells) 


r c 


12.5 um 


Position of Casparian strip 


r e 


40 um 


Root radius 


To 


500 uM um" 2 s" 1 


Maximal amount of ZIP 


Ho 


5 uM um -2 s _1 


Amount of HMA4 


i 


0.01 - 1 


Scaling factor for regulatory time scale 



Table 2: Cell sizes used in simulations. Values 
correspond to a typical young Arabidopsis thaliana 
root. 

Cell type Left radius Right radius 

(H (um) 

Pericylce 6 9.75 

Endodermis 10.25 14.75 

Cortex 15.25 29.75 

Epidermis 30.25 40 
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FIGURES 



Figure 1: Scheme of root tissues. ZIP transporters are localized on the epidermis (ep), the cortex 
(co) and the endodermis (en). Endodermis cells are surrounded by the Casparian strip (cs). HMA4 
transporters are localized in the pericycle (pc). The symbols r e , r c , and r x denote the radial 
coordinates at the outisde of the epidermis, the casparian strip and the inside of the pericycle, 
respectively. Arrows illustrate the direction of transport. 

Figure 2: Flow velocity and volume fraction. Relative flow velocity (top) and volume fraction 
(bottom) as functions of the radial coordinate. Tissues include the epidermis (ep), cortex (co), 
endodermis (en) and pericycle (pc). 

Figure 3: Scheme of ZIP regulation model. The activator^ dimerizes and induces gene activity, 
mRNA transcription and translation of the transporter protein. Zinc is transported across the cell 
membrane and causes the inhibitor I to inhibit^. Figure adapted from Claus and Chavarria-Krauser 
(2012). 

Figure 4: Spatial distribution of zinc in wild type. Symplastic (solid line) and apoplastic (dotted 
line) zinc concentration for high (25 uM; top) and low (1 uM; bottom) medium concentration. In 
addition, grey boxes illustrate the ZIP activity in epidermis (ep), cortex (co) and endodermis (en). 

Figure 5: Variation of HMA4 and flow velocity. Symplastic zinc concentrations are shown along 
the radial coordinate in epidermis (ep), cortex (co), endodermis (en) and pericycle (pc) cells. Top: 
HMA4 levels were increased to two (dashed line) or three times (dotted line) compared to the wild 
type (WT) level (solid line), while the flow velocity q was kept constant at 4 um s" 1 . Bottom: Flow 
velocity q was decreased from 4 um s" 1 (solid line) to 2 um s" 1 (dashed line) and 1 um s" 1 (dotted 
line), while HMA4 was kept constant at the WT level. 
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Figure 6: Symplastic zinc gradient. Average symplastic zinc gradient shown as a function of 
HMA4 in multiples of the wild type levels (solid line) and of the external zinc concentration Z e 
(dashed line). Negative numbers reflect a gradient oriented towards the stele. 

Figure 7: Change in external zinc concentration. Starting from the steady state distribution of 
zinc in the symplast (solid line) and apoplast (dashed line) at an external concentration of 1 uM, the 
medium concentration was changed to a high zinc condition (10 uM). Grey boxes show ZIP activity 
for epidermis (ep), cortex (co) and endodermis (en) cells. The new steady state is reached at t = 20 
s. 

Figure 8: Change in flow velocity. Starting from the steady state distribution of zinc in the 
symplast (solid line) and apoplast (dashed line) at q = 0.05 |j,m s" 1 , the flow velocity was increased 
to qo = 4 urn s" 1 . Grey boxes show ZIP activity for epidermis (ep), cortex (co) and endodermis (en). 
A new steady state is reached after 5 s. External zinc was kept constant at Z e = 25 uM. 

Figure 9: Variation of regulation time scale. The time scaling factor x of the ZIP regulation model 
was decreased from x = 1 (red line) to x = 0.1 (blue line) and x = 0.01 (green line) and time courses 
were simulated for a change of the external zinc concentration from 1 |aM to 10 uM. A, the 
transition between the two steady states (black squares) is shown as a phase diagramm of ZIP 
activity again symplastic zinc. For x = 1, the system reaches the new steady state after a minor 
overshoot. For x = 0.1 and x = 0.01, the new steady is unstable and concentrations oscillate on a 
limit cycle. B, time course of the symplastic concentration showing clearly the oscillations. 
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Summary 

A cell-based model for the uptake and transport of zinc in roots of Arabidopsis 
thaliana and Arabidopsis halleri is derived. The model consists of a coupled sys- 
tem of ordinary differential equations describing the regulation of ZIP transporters 
and 1-D partial differential equations describing the transport in the symplast and 
apoplast. It considers, thus, the internal structure of the root tissue and couples 
transport phenomena with regulation networks. A system of ordinary differential 
equations for the xylem is also derived from a transport model and coupled via a 
boundary condition to the 1-D model of the apoplast. 

S.l Overview 

Transport in the root is assumed to be mostly in radial direction. This allows to reduce 
the three dimensional problem into coupled one dimensional radially oriented problems. 
The structure along the radius is shown schematically in Fig. 1 of the main manuscript. 
The root is assumed to be composed of following cell types (from outside to inside): 
epidermis (ep), cortex (co), endodermis (en) and pericycle (pc). The cell layers extend 
from radius r x to r e . Surrounding the endodermis cells, a perfectly unpermeable Cas- 
parian strip (cs) at position r c is assumed. Epidermis, cortex and endodermis cells are 
allowed to have ZIP transporters (influx), while pericycle cells have only HMA4 efflux 

* andres.chavarria@bioquant.uni-heidelberg.de 
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Figure S.l: Scheme of the domains correspoding to apoplast fii, ^2 and $^4, and symplast 
The root is assumed to be composed of a periodic assembly of such 
domains. 

transporters. Following the results of Talke et al. [2006], the expression of HMA4 is as- 
sumed to be independent of the zinc concentration and will be included into the model 
as a given efflux. Both transporters ZIP and HMA4 are assumed to be saturable and to 
follow Michaelis-Menten kinetics. The expression of ZIP in the epidermis, cortex and 
endodermis is allowed to adapt to the current internal zinc status based on the dimerising 
activator-inhibitor model proposed by Claus and Chavarria-Krauser [2012]. Depending 
on the average internal zinc concentration, each cell decides independently of the others 
its expression level of ZIP. The resulting amount of transporters is assumed to be evenly 
distributed on the plasma-membrane and will vary in general from cell to cell (but not 
within). This assumption is supported by the HMA2 expression pattern found in A. 
thaliana, Sinclair et al. [2007] (Fig. 1(c) of that publication). Note here that, the model 
uses no other type of signal beside the internal zinc concentration. Hence, coordination 
is achieved merely by zinc fluxes. 

To reduce the problem to a system of 1-D partial differential equations, some as- 
sumptions on the geometry are needed. First, the root is best described by cylindrical 
coordinates and reduction will be achieved by averaging over the axial coordinate z and 
the azimuth <p. Second, assume that the layer of interest is made up by periodic repeti- 
tion of Cl := f^i U O2 U Q3 U in azimuthal and axial direction. Here, Q.3 denotes the 
symplast and Oi, Q2, and O4 different parts of the apoplast respectively and as indicated 
in Fig. S.l. And third, assume that the contribution of the apoplastic "edges" O4 to the 
overall zinc uptake is negligible. This assumption is sustainable, as O4 is thin compared 



2 



to the rest. Each cell layer presented in Fig. 1 could be described by such a periodic 
assembly. We will set all Qi to span over all cell layers from r x to r e . 

Cells have a complex internal structure with organelles, such as vacoules, nucleus, etc. 
They are also interconnected by plasmodesmata - channels which traverse the cell wall 
and represent a substantial reduction in flow cross section. To avoid the treatment of 
these internal structures, we regard the cell content to be a sort of porous medium with 
a given volume fraction, i.e. the vacoules are not treated as separate structures. Cell 
walls are also best described as a porous medium of constant structure and porosity. 
Therefore, we introduce a volume fraction = $(r) for the symplast, which depends 
only on the radial position. This assumption is valid in view of the periodic structure of 
the cell and orientation of cell layers (Fig. 1). The volume fraction of the apoplast will 
be denoted by the constant 

S.2 Regulation 

Regulation is modelled as in Claus and Chavarria-Krauser [2012] by an ordinary differ- 
ential equation system. One system per each cell type listed in Fig. 1 is needed 

1G (K A 2 a (l - G a ) - G a ) , 
7M (G a - M a ) , 

7t (M a — T a ) , for t G (0, oo), a = en, co, ep , (SI) 

-f A (l-TA a I a -A a ) , 

11 {Tl (a ~ r' A a I a -(1 + T 7 C Q ) la) , 

where G a is the gene expression level, M a the transcript level, T a the transporter level, 
A a an activator and I a an inhibitor, and £ a the internal zinc concentration. The 7 
factors are related to the time scales of each single reaction step. This system was 
non-dimensionalized in such a way that all variables assume values between and 1. 
Pericycle cells are not included, because these efflux cells can be assumed to not express 
an influx transporter. Fig. 3 in the main manuscript presents a scheme of the processes 
modelled by Eq. (SI). 

The model expects one single value for the internal zinc concentration ( a (t), while 
the transport model delivers a distribution Z = Z(x,t) with x G O3. Thus, the zinc 
concentration varies inside a single cell. To circumvent this issue, we assume that a cell 
senses the average zinc concentration inside it 

Ca(t) = (q 1 ^ c ^ Z(x, t) dx for t G [0,oo), a = en, co, ep , (S2) 



~dT 
dM a 

dt 
dT a 

dt 
dA a 

dt 
dl a 
dt 
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where is a scaling factor (£ a is non-dimensionalized) , integration is over the cell 
C a C ^3 and n(C a ) is its volume. The transporters regulated by Eq. (SI) are assumed 
to be distributed evenly on the surface dC a of the cell, and a distribution of transporters 
is constructed as follows 

ep 

T(x,t) = T *(t)Xdc a (x) , fOT (M) e H 3 x [0,00) (S3) 

ct=en 

with the characteristic function 

1 if x € <9C a , 
if x g 3C a . 

HMA4 efflux transporters at the pericycle are included in a similar manner 

H{x) = H pcXd c P M , xgH 3 , (S4) 

where the assumption that the level H pc of HMA4 is constant was used, Talke et al. 
[2006]. 

S.3 Water flow 

Before treating our main topic of zinc transport in the root, we construct a simple model 
for water flow. Membranes restrict the movement of zinc, but water carries it along the 
flow path (advection). This process influences the distribution of zinc and determines 
how fast variations in external zinc concentration spread in the system. Ultimatelly, 
advection is essential to the regulation patterns. 

To avoid a complete treatment of water fluxes in root tissues, we focus only on mass 
conservation delivering the flow speed by consideration of effective flow cross sections. 
Variation of cross section in the symplast is included through the volume fraction <J>(r). 
Water fluxes in the apoplast are assumed to be small compared to those in the symplast 
and will be neglected here. This includes the assumtpion that there is no exchange 
between symplast and apoplast, although the apoplast is believed to contribute to the 
total flux, Steudle [2000]. We assume that epidermal cells take up water from the medium 
with a constant flux density go- During its pathway to the xylem, water is assumed to be 
conserved and to have a flow speed that depends only on volume fraction and geometry 
(radial convergence). This approach is very simplistic and other more sophisticated 
models have been proposed, for example Katou and Furumoto [1986], Katou et al. [1987], 
Taura et al. [1988], Katou and Taura [1989], Murphy [2000]. Modeling water fluxes in 
plant tissues is a complex problem which deserves treatment of its own and is out of 
scope of this manuscript. 

Mass conservation for an incompressible fluid reads in the symplast 

div($v) = in Q 3 , (S5) 



Xdc a ( x ) = 
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where $ is the volume fraction and v is the flow velocity. Eq. (S5) is expressed in 
cylindrical coordinates (r, ip, z) to reflect the geometry of roots 

^d r (r^v r ) + ^-d^v ip )+d z (^v z ) = . (S6) 

where v r , v v , and v z denote the radial, azimuthal, and axial component of the velocity 
v. To the end of reducing the model to 1-D we consider the surfaces 

r,(r) :={xeOi|^ + ^= r 2 } for i = 1, 2, 3 , (S7) 

which can be described in cylindrical coordinates by 

(r, <p, z) £ {r} x (0, ^ ,i) x (0, z ,i) for i = 1, 2, 3 , 

with azimuth ipQi and height 2o,i of the considered domain. Note the use of different 
polar coordinate systems for Q2, and Let ^ (r) = r <£>o,i £o,i denote the area of 
Ti (r) . Introduce the averaged radial velocity in the symplast 

v(r) := — l — [ v r (r,Lp,z)d-f , for r G [r x ,r e ] , 
M3 (,?•) Jr 3 (r) 

and consider the corresponding average of Eq. (S6) over T3 (r): 

— l — [ -d r (r<S>v r )d 7 = l — f (-O^v^ + d^v^dj , 

Mr) Jr 3 (r) r Mr) J rs{r) \r J 

where the second and third terms in Eq. (S6) were moved to the right hand side. With 
the surface element d7 = r dtp dz the left hand side of the equation is 



f 1 1 

/ -dJr <S>v r ) d7 = -dJr <&v) , 

Mm r r 



M r ) Jr 3 (r) r 

while the terms on the right hand side correspond to a two dimensional divergence and 
can be transformed into an integral over the boundary dT^r). This boundary integral 
is zero, based on the assumption that the apoplast and symplast do not exchange water 
(consequence of assuming no water fluxes in the apoplast). We obtain an equation for 
the average flow velocity in the symplast 

d r (r $ v) = for r x < r < r e , (S8) 

Hr=r e = »> ( S9 ) 

where the water influx qo was prescribed on the root surface (r = r e ). This system can 
be solved by integration rendering 



5 



Fig. 2 in the main manuscript shows the flow velocity used in the simulation. The 
volume fraction in the symplast <J>(r) considers the periodic space restriction produced by 
vacoules and plasmodesmata. An estimation based on Rutschow et al. [2011] delivers a 
volume fraction in plasmodesmata of the order of 0.15, while the vacoule was assumed to 
make up 80% of the cell volume ($3 = 0.2). Table S.2 presents the geometry parameters 
on which the calculations were based. This parameter set represents a typical A. thaliana 
root. 



S.4 Zinc transport 

Having determined the flow velocities in the compartments considered, we are able to 
move on to the task of finding a model for zinc. Its transport can be modelled by a 
diffusion advection problem, that states the conservation of zinc 

dt{^Zi) - dw(^D gradZi) = in ^ x (0, 00) , i = 1, 2 , (Slla) 
dt($Z 3 ) + div (&vZ 3 -$D grad-Z 3 ) = inft 3 x(0,oo), (Sllb) 

where ^ is the volume fraction of the apoplasts Qi and Q2, Zi is the zinc concentration 
and D a diffusion coefficient. Transport in the apoplast is assumed to take place only 
by means of diffusion, as the water flow velocity was assumed to be zero. 



S.4.1 Reduction to 1-D 

Define Z{ as the average of Zi over Tj (r): 



Zi(r,t) := —— [ Zi(x,t)dj , for (r, t) G [r x , r e ] x [0, 00) . 

We proceed to average Eqs. (Slla) and (Sllb). We will average exemplarily Eq. (Sllb). 
The result will apply also to Eq. (Slla) by exchanging <E> with and setting v = 0. 
Treatment of the time derivative term is straightforward 



^— f d t ($Z 3 ) d 1 = d t {<S>Z 3 ) 



m r ) Jr,(r) 

The contribution of the advection term in Eq. (Sllb) is 

— l — f div($vZ 3 )d 7 = 1 d r (r$— l — f v r Z 3 d 1 )+-^[ $«„Z 3 ds, 
H3{r) Jr 3 (r) r y /i 3 (r) 7r 3 (r) / MsW J d r 3 ( r ) 

where v n is the normal velocity on dT%{r). For the sake of generality the integral over 
dT 3 {r) is kept, although it is zero by the assumption of no water exchange between the 
apoplast and symplast (v n = 0). By the mean value theorem, there exists a v (r) so that 



/ v r Z 3 dj = v (r) / Z 3 d-y . 
Jr s (r) Jr 3 (r) 
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This renders the approximation 

— l — [ div($vZ 3 )d 7 » -d r (r<f>vZ 3 ) + — l — [ $v n Z 3 ds , 
mr) Jr 3 (r) r n 3 (r) J d r 3 ( r ) 

where v ~ v was assumed. In the last step we had to make a somewhat crude approx- 
imation. However, without explicitely solving the whole 3-D problem, this is the best 
that can be done. The contribution of diffusion is obtained by expressing the divergence 
and gradient operators in cylindrical coordinates. The terms containing an r-derivative 
are 

^- / -d r {<S>Drd r Z 3 )d 1 = --d r (<S>Drd r Z 3 ) . 

The term containing ip- and z-derivatives are transformed into an integral over the 
boundary of T 3 (r) 

[ (-d v (—d v Z 3 )+d x ^Dd x Z 3 ))dry = ^- / $Dd n Z s ds, 

Jr 3 (r) V V r J J to(r) J 9 r 3 (r) 



where d n Z 3 is the normal derivative of Z 3 on dF 3 (r). Summation of this boundary 
integral with the one obtained from the advection term delivers the average flux through 
the boundary 

— l — [ $(v n Z 3 -D d n Z 3 ) ds = -Q 3 , 

2 (S12) 
Qi ■= ^2^ij(r)Jj(r) 

3=1 

where the flux densities Jj(r) represent the transport through ZIP and HMA4 (developed 
in next section), Q 3 denominates the source of zinc, and Oij is 

/ -2/20,1 \ 
K)(r) := -2/r^ . (S13) 

\ 2/zo,3 2/r99 ,3 / 

The sign convention of J, L is as follows: positive for a flux from the apoplast into the 
symplast and negative vice versa, crjj is related to the ratio of the length of the pieces 
composing the boundary dT 3 (r) to the area of T 3 (r) and weights, thus, the flux densities 
Jj to account correctly the change in average concentration. 

Putting everything together, the following reduced model for the symplast is ob- 
tained 

d t (<S>Z 3 ) + ^-d r (r<S>vZ 3 -<S>Drd r Z 3 )=Q 3 in (r x , r e ) x (0, 00) . (S14) 
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It is practical for the implementation of a numerical algorithm to rewrite the equation 
into a conservative form by introducing rZ; as a variable. Applying the results for 
to fli and Q2, the following system is obtained 

dt(^frZ 1 )+dr((D/r)9rZ 1 -'9Ddr{rZ 1 )) =rQ 1 , 

d t ^rZ 2 ) + d r ((D/r)^rZ 2 -^Dd r (rZ 2 )) = rQ 2 , in (r x , r e ) x (0, oo) . (S15) 
d t ($rZ 3 )+d r ((v + D/r) <P rZ 3 - <t> D d r (rZ 3 )) =rQ 3 , 

Initial values and boundary conditions are developed in Section S.4.3. 
S.4.2 Flux densities Jj 

The zinc sources Qi on the right hand side of Eq. (S15) depend on the flux densities Jj, 
which still need to be specified. The zinc flux through transporters can be modelled by 
a saturable pointwise reaction mechanism 

Zi + T Z-T — >Z 3 + T, fori = 1,2, 
Z 3 + H Z-H ^Zi + H . 

We will assume that the above reactions follow Michaelis-Menten kinetics and introduce 
as in Claus and Chavarria-Krauser [2012] a saturation function 

f{z ' K) = zTK' (S16) 

where K is the corresponding Michaelis-Menten constant. We introduce versions of 
T(x,t) and %(x) which depend on the radius by averaging over dT 3 {r) 

T ( r ' *) := /J ? ^^ / *) ds > for ( r > *) G Vx, re] x [0, oo) , (S17a) 

H{dr 3 (r)) Jar 3 (r) 

Note that, T(r,t) is equal to T en (t), T co (t) and T ep (t), respectively, for an r inside one 
of these cells, zero elsewhere. H(r) is equal to H pc for an r inside the pericycle and zero 
elsewhere. The reaction probability depends on ^>T(r,t) instead of only T{r,t), because 
Qi and are porous media and only the reduced amount ^/T(r, t) has actually contact 
to Z\ and Zi- No correction is needed for H(r), as the cytoplasm can be assumed to 
have direct contact with the membrane, so that the complete H(r) can react with Z 3 . 
In total, the flux densities Jj are modelled as 

Jj{r,t) = $T T{r,t)f{Z j {r,t),K t )-H Q H(r)f(Z 3 {r,t),K h ) , j = 1,2 , (S18) 

where To and Hq are constants that characterize the true amount of transporters (non- 
dimensionalized regulation) . 
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S.4.3 Initial values and boundary conditions 



Eq. (S15) needs initial values and suitable boundary contions to obtain a well posed 
(solvable) problem. The apoplast is assumed to have access to a perfectly stirred medium 
of concentration Z e . A concentration Z x (t) is prescribed at the xylem. This concen- 
tration depends on the flux of zinc through HMA4 and a model will be developed in 
the next section. The impermeability of the Casparian strip is considered by setting a 
no-flux condition. In total, we prescribe for the apoplasts Q\ and Q2 



(D/r)^rZi - ^Dd r (rZi 
(D/r)^rZi -VDd r {rZi 



z x (t) 
, 
, 
z e , 



for t £ [0,oo), 2 = 1,2 



(S19) 



The zinc flux is prescribed at the boundary of the symplast O3 

(v + D/r) <Z>rZ 3 -<S>Dd r (rZ 3 ) = r H H pc f (Z 3 (r x ), K h ) 

r=r x 

(v + D/r) <S>rZ 3 -$Dd r (rZ 3 ) = -r <S>T T ep f (Z e , K l ) , 

r=r e 

for t G [0, 00) . 



(S20a) 
(S20b) 



Published experimental results that capture the dynamics of regulation, focus on 
changes from one steady state at a given external concentration to another steady state 
for a different concentration (e.g. zinc resupply). Therefore, the initial conditions used 
here are solutions of the stationary version of Eq. (S15) 



Zi(r,t) 



t=o 



Zi(r) for r G [r x ,r e ] , i = 1,2,3, 



(S21) 



where the Z% fulfill one of the following equations 



d r ((D/r)VrZ i -yDd r (rZ i ))=rQ i , in (r x , r e ) x (0, 00) , % = 1, 2 , (S22) 
d r ((v + D/r) $rZ 3 -$Dd r (rZ 3 )) = r Q 3 in (r x , r e ) x (0, 00) , (S23) 

with boundary conditions Eqs. (S19), (S20a) and (S20b). In general, the steady state for 
given Z e and qo will be determined, used as a initial condition, one of these parameters 
changed, and the dynamics of the transition captured. 



S.4.4 Xylem 

A calculation of the apoplastic zinc concentration in the region enclosed by the casparian 
strip (r x < r < r c ) needs the concentration of zinc in the xylem. For simplicity, we will 
pose a model for the central cylinder < r < r x (i.e. stele without the pericycle) 
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and account the true size of the xylem by a constant volume fraction *f> x . The domain 
describing this tissue will be denoted as Q x , where the x stands for xylem. 
Eqs. (S5) and (SI lb) apply also to this tissue 

div(tf x v x ) = in n x , (S24) 

d t (y x Z x ) + div(y x v x Z x -y x D grad^) = in Q x x (0, oo) , (S25) 

where v x is the flow velocity and Z x is the zinc concentration in the central cylinder £l x . 
An average can be obtained as in Sees. S.3 and S.4. The main difference is that the 
surface over which the average is created is here 

T x (z) := {x£n x \x 3 = z} . 

Conservation of water delivers 

d z (* x v x ) = -—*{r x )v{r x ), (S26) 

fx 

where ®(r x )v(r x ) is the flow velocity of the water being delivered from the symplast, 
and the average velocity in the xylem is defined as 

v x (z) := , . / v Z:X (r,ip,z)dj , forzG[0,L] , 

where L is the length of the root portion considered. Expression of &(r x ) v(r x ) by Eq. 
(S10), integration of Eq. (S26), and assumption of v x (0) = delivers 

Vx (z) = -^ r -±^- z for*€[0,L], (S27) 

r x fx ^ X 

which is a linear function of z. Remeber that go < so that v x (z) > for z > 0. Eq. 
(S27) is based on the assumption that qo is constant, which will not be true in reality. 
The pressure gradient between the xylem and the medium will fall with z and, hence, 
v x (z) cannot grow linearly indefmitelly and will stagnate at a constant value. However, 
v x (z) will behave similar to Eq. (S27) in a region near z = 0. We focus on this region 
and assume validity of Eq. (S27). 

An average for the equation describing the conservation of zinc is obtained readily 

dt(* x Z x ) + d z (* x v x Z x -* x Dd z Z x ) = -H H pc f(Z 3 (r x ,t),K h ) , (S28) 

r x 

where the boundary condition Eq. (S20a) divided by r was used and the average zinc 
concentration is defined as 

Z x (z,t) := 1 / Z x (x,t)d 7 for (z, t) G [0, L] x [0, oo) . 
M(r*(«)) Jr x (z) 
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A preliminary simulation of this equation with no-flux and open vessel conditions at 
z = and z = L, respectively, shows that Z x is almost constant in space. Hence, we set 
Z x (z,t) Z x (t), use that ^> x is constant and express v x by Eq. (S27) to obtain 



dZ x 1 2 



(^q Z x + H H pc f(Z 3 (r x ,t),K h fj forte (0,oo) , 



dt * x r x \r x ™ u V ~" (S29) 

z lt=o _ Z o • 

Note again that qo < 0, so that this equation has a non-trivial positive steady state 
solution 

T = - TxI ^f(Z 3 (r x ),K h ) . (S30) 
r e qo 
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Figure S.2: Scheme of the steps involved in coupling the numerical schemes. 
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A 2 v = A ^ o 



Figure S.3: Scheme of effect of cylindrical geometry on water flow velocity. The flow 
cross section becomes smaller for smaller radii resulting in acceleration of 
the flow: v = ^ go > qo for A\ > A2. 



Table S.l: Parameters used in the regulation model. Values taken 
from Claus and Chavarria-Krauser [2012] 
Parameter K l * K T V Tj jc 7m It 1A 
Value 13 fiM 20 38 167.2 1000 4 4 1 1 

* Value for ZIP1, Grotz et al. [1998] 
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Table S.2: Geometry parameters used in the simulation. These cor- 
respond to a typical Arabidopsis thaliana root. 
Parameter r x r c r e z ,i ^0,3 <Po,3 

Value 6 fim 12.5 \xm 40 fim 0.5 fim 135 /im 7r/10 



Table S.3: Further parameters used in the simulation. 
Parameter Value Description 

K h 1 fj,M Michaelis-Menten constant for HMA4 

Co 166.67/uM Scaling factor that dimensionalies the internal zinc concentration 
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